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Abstract 

Dendritic cells (DCs) are critical for regulating CD4 and CDS T cell immunity, controlling Thl, Th2, and Thl7 commitment, 
generating inducible Tregs, and mediating tolerance. It is believed that distinct DC subsets have evolved to control these 
different immune outcomes. However, how DC subsets mount different responses to inflammatory and/or tolerogenic 
signals in order to accomplish their divergent functions remains unclear. Lipopolysaccharide (LPS) provides an excellent 
model for investigating responses in closely related splenic DC subsets, as all subsets express the LPS receptor TLR4 and 
respond to LPS in vitro. However, previous studies of the LPS-induced DC transcriptome have been performed only on 
mixed DC populations. Moreover, comparisons of the in vivo response of two closely related DC subsets to LPS stimulation 
have not been reported in the literature to date. We compared the transcriptomes of murine splenic CDS and CDllb DC 
subsets after in vivo LPS stimulation, using RNA-Seq and systems biology approaches. We identified subset-specific gene 
signatures, which included multiple functional immune mediators unique to each subset. To explain the observed subset- 
specific differences, we used a network analysis approach. While both DC subsets used a conserved set of transcription 
factors and major signalling pathways, the subsets showed differential regulation of sets of genes that 'fine-tune' the 
network Hubs expressed in common. We propose a model in which signalling through common pathway components is 
'fine-tuned' by transcriptional control of subset-specific modulators, thus allowing for distinct functional outcomes in closely 
related DC subsets. We extend this analysis to comparable datasets from the literature and confirm that our model can 
account for cell subset-specific responses to LPS stimulation in multiple subpopulations in mouse and man. 
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introduction 

Dendritic cells (DCs) are key regulators of T cell responses. DCs 
are essential for priming naive T cells and are also believed to 
control their effector fate. The DC lineage can be subdivided into 
multiple distinct subsets, some of which show intrinsic functional 
differences that are known to drive distinct immune outcomes [1]. 
However many subset-specific functional differences remain 
poorly understood. Here we have used a global systems approach 
to DC function as a means of exploring their distinct in vivo roles 
in the immune response. 

Advances in systems biology have clearly demonstrated that 
linear signalling cascades poorly represent the complexity of 
immune signalling (reviewed in [2]). Rather than comprising linear 
pathways, immune signalling involves interactions between 
thousands of distinct proteins communicating within a complex 
network. These networks are organised by a set of highly 



connected proteins (known as Hubs) that are essential for receiving 
and distributing multiple signals within the network [3-5] . Due to 
their key role in the connectivity of complex signalling networks. 
Hubs both reflect mechanism and provide biomarkers for cell 
types and signalling events [3-5]. It is not yet known whether 
differences in Hub usage contribute to cell-specific differences in 
signalling networks. 

In vivo toll-like receptor 4 (TLR4)-dependent responses to 
bacterial lipopolysaccharide (LPS) provide an ideal model in which 
to test whether closely related cell subsets show differences in their 
immune signalling networks, since a wide array of cell types 
express TLR4 and respond to LPS [6-9]. Studies using systems 
biology approaches to investigate LPS responses have primarily 
focused on clarifying shared mechanisms rather than defining the 
differences between closely related cell subsets. Published studies 
have shown that LPS responses are initially propagated through 
two sets of adaptor molecules: ToU-interleukin- 1 receptor (TIR) 
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domain-containing adaptor protein (Tirap) and Myeloid differen- 
tiation primary response 88 (Myd88) {Tirap-Myd88, the 'Myd88- 
dependent patliway'), or TIR-domain-containing adapter inter- 
feron-P-inducing factor (Trif) and TRIF-related adaptor molecule 
(Tram) (Trif-Tram, the 'Myd88-independent pathway'). Addition- 
ally, a set of Hubs responsible for orchestrating signaUing 
outcomes in response to LPS has been defined [6-9] . These Hubs 
are essential for signal propagation and belong primarily to the 
tumor necrosis factor receptor associated factor (TRAF), interleu- 
kin-1 receptor- associated kinase (IRAK), mitogen-activated pro- 
tein kinase (MAPK) and nuclear factor kappa-liglit-chain-enhanc- 
er of activated B cells (NFkB) families (extensively reviewed in [6- 
9]). 

Transcriptional analysis of the LPS response in murine DCs has 
generally been confined to cells differentiated in vitro from bone 
marrow precursors, with a single report of the response of 
unfractionated ex vivo splenic DCs [10-15]. The splenic DC 
compartment comprises distinct cell subsets expressing either CD8 
or GDI lb and manifesting different basal transcriptional pro- 
grams [16-26]. CD8 DCs are thought to uniquely cross-present 
antigen to CD8 T cells and are the major producers of interleukin 
(IL-12) for the regulation of Thl responses, while GDI lb DCs are 
thought to be dominant in the regulation of CD4 T cell responses 
and Th2 immunity [1,24,27], although not all models support 
these functional distinctions [28]. Inactivation of key transcription 
factors, including IRF8, BATF3, IRF4 and Ikaros, selectively 
interferes with development of CD8 or GDI lb DC subsets 
[29,30]. However differences in inflammatory signalling pathways 
in the two subsets remain poorly defined. Both have been reported 
to respond direcdy to LPS in vitro and in vivo, although they 
express relatively low basal levels of TLR4 [31,32]. In vitro 
stimulation with LPS induces equivalent production of tumor 
necrosis factor alpha (TNFot) and IL-6, but higher production of 
IL-12 in CD8 DCs, suggesting that while both subsets share 
common signalling pathways and TLR4 potency, subset-specific 
differences are also present [31]. Both DC subsets also respond to 
a number of mediators released by DCs and other cell types in 
response to LPS, so that their response to in vivo LPS 
administration comprises a network of direct and indirect effects 
that joindy control their ability to differentially stimulate T cells 
[1,27]. Such effects may not be adequately modelled by in vitro 
LPS stimulation of purified cells. 

In this study, we administered LPS in vivo and isolated splenic 
DCs from untreated and LPS-treated mice with minimal 
manipulation. We then used RNA-Seq of flow-sorted samples 
pooled from multiple mice to compare the physiological responses 
of closely related DC subsets to in vivo LPS exposure. Using a 
hyper-stringent method for choosing differentially expressed genes 
in the DESeq R package [33], we show that CD8 and GDI lb 
splenic DC subsets respond differendy to in vivo LPS stimulation, 
and that many of the transcriptional changes previously defined in 
the LPS response of unfractionated DCs are present in only one of 
the two subsets. 

We used network analysis to identify a subnetwork in each DC 
type in order to elucidate the mechanisms underlying the observed 
differences. Such subnetworks are thought to contain key 

regulators of the measured response. Importantly, they are self- 
reinforcing, and are therefore less susceptible to variability in 
individual gene detection [3-5]. The 2 DC subsets generally 
expressed the same set of core LPS response molecules, many of 
which served as Hubs in each subset-specific subnetwork. Both 
subsets also expressed a common set of cell-surface receptors 
required for responses to secondary mediators released after LPS 
stimulation. However, the sets of proteins interacting with these 



core Hubs were significantiy different in the 2 subsets. Important- 
ly, the majority of such interacting proteins, including Atf3, 
TnfaipS (A20), Tradd and Cdknla, are already known to be 
modulators of common signalling pathways, although they had not 
previously been accorded subset-specific roles. These data support 
a model in which distinct immune responses to the same stimulus 
are achieved by differential 'fine-tuning' of core pathways by 
subset-specific modulators. Finally, we validated our hypothesised 
model using meta-analyses of other cell populations, showing its 
relevance to inflammatory LPS signalling in multiple cell subsets in 
mouse and human. 

Results 

Differential activation of DC subsets by in vivo exposure 
to LPS 

Spleen cells were harvested from steady-state (n = 5) and LPS- 
treated (n = 10, 24 hours after 25 \ig LPS i.p) mice and DC subsets 
purified using magnetic bead enrichment for lineage (CD 19, B220, 
CD3, Gr-1, Terl 19)-negative, CD 1 1 c-positive cells, followed by 
flow sorting according to the gating strategy shown in Figure lA. 
We confirmed that DCs were indeed activated by in vivo LPS 
administration by comparing expression of the activation marker 
CD86 to that of steady-state cells (Figure 2A-B). As expected [34], 
both DC subsets responded to LPS in vivo by up-regulating CD86. 
We also confirmed that our enrichment and gating strategy for 
steady-state and LPS-treated DCs excluded monocyte-derived 
DCs identified on the basis of coexpression of FcyRl (CD64) and 
FcsRlot [35] (Figure SI). RNA-Seq was then performed on the 4 
RNA samples, using standard techniques as described in the 
Materials and Methods section. Significantly differentially ex- 
pressed genes were identified using a hyper-stringent method from 
the R package DESeq at a p-value cut-off <0. 05 [33]. We assessed 
the quality of our RNA-Seq data by comparing our steady-state 
data with pubhshed results from the literature. We first compiled a 
set of the top 50 prototypical subset-specific genes, based on 
published studies of mRNA and protein expression, and showed 
that our steady-state data faithfully recapitulated the published 
patterns (Figure IB). Next we calculated the overlap between our 
steady-state data and 9 published datasets (datasets 1-9 listed in 
Table SI [16-26]) using a hypergeometric test (Figure IC, see 
Materials and Methods). The highly significant overlap between 
differentially expressed genes in our steady-state dataset and each 
of the published datasets indicated that our steady-state data were 
remarkably consistent with previously published data. These 
analyses indicated that our dataset provided a suitable measure 
of gene expression, even though it contained only a single pooled 
sample for each experimental condition. 

We next analysed LPS-induced gene expression by comparing 
LPS stimulated with steady-state data for each subset. LPS 
stimulation of CD8 DCs led to a significant change in the 
expression of 481 genes (397 upregulated, 84 downregulated), 
while in CD lib DCs there was a significant change in the 
expression of 471 genes (428 upregulated, 43 downregulated) 
(Table S3). We also confirmed that both subsets expressed 
detectable Tlr4 mRNA (CD8 DCs: 59 and CD lib DCs: 92 
normalised counts per million). Interestingly, the observed LPS 
responses were highly subset-specific. Thus 49'/o of differentially 
regulated transcripts in CD8 DCs were subset-specific, while the 
corresponding figure for CDl lb DCs was 48% (Figure 2C). Many 
of these differentially regulated genes are known to be important 
immune effector genes, suggesting an LPS-regulated functional 
divergence between the 2 DC subsets. Compared with CD8 DCs, 
CDl lb DCs selectively upregulated a wider range of transcripts 
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Figure 1 . Comparison of steady-state spleen DC subsets. (A) Steady-state splenic DCs were magnetic bead-enriched for CD1 9"B220"CD3"Gr- 
I^Terl 19"CD11c* cells. MHClTCDI 1c* cells (circled) were then sorted for CDS (blue gate) and CD1 lb (orange gate) subsets and RNA prepared and 
analysed by RNA-Seq. (B) Heatmap showing the relative expression of the 50 most commonly defined and validated markers for CDS and CDIIb 
subsets. Data are presented as fold changes (CD1 Ib/CDS), all of which were statistically significant with an associated p-value <0.05. Orange denotes 
genes that were increased in CD1 1 b DCs while blue denotes genes that were decreased in CD1 1 b DCs (and thus increased in CDS DCs. (C) Overlap 
between genes that were significantly differentially expressed between CDS and CD1 1 b DCs in our dataset and in 9 previously published microarray 
datasets derived from splenic DC subsets (datasets 1-9 listed in Table SI). The significance of overlap between the gene list from our dataset and 
those from each of the published datasets was calculated using a hypergeometric test to assess the consistency/quality of our results. Data are 
presented as 1 /p-value on a log scale with all overlaps reaching a significance cut-off <0.05. 
doi:1 0.1 371 /journal.pone.01 0061 3.g001 



encoding cytokines/chemokines (including Ccl6, Ccl7, Ccl22, 
Cxcl3, CxcllO, Cxclie, Illb, 1115, and Tnf), while CDS DCs 
selectively upregulated co-stimulatory molecules such as Cd274 
(PD-Ll), Icosl, and Tnfsf4 (OX40L) (Figure 2D). These transcrip- 
tional changes are likely to be mediated by a combination of direct 
(LPS-TLR4) signals and secondary soluble (cytokines/ chemokines) 
and/or cell-to-cell signals from the LPS-activated splenic micro- 
environment. 

To identify known biological pathways underlying the differ- 
ences between the subsets, we performed gene ontology (GO) term 
over-representation analysis on the significantly differentially 
regulated genes. The GO term analysis afforded further evidence 
that our data from single pooled samples provided a valid measure 
of gene expression. Thus genes annotated by GO terms associated 
with LPS-stimulation (response to lipopolysaccharide and cellular 
response to lipopolysaccharide) and with the general inflammatory 
response (cytokine-mediated signalling pathway, immune re- 
sponse, inflammatory response, innate immune response) were 
significantly enriched in both subsets, although the individual 
genes within the GO categories differed between the subsets 
(Table 1). GO terms associated with regulation of the apoptotic 
process were enriched in both subsets although they reached 
statistical significance only in GDI lb DCs, while the GO term 



'negative regulation of the inflammatory response' was signifi- 
cantly enriched only within the CDS DC subset (Table 1). 
Interestingly both apoptosis and negative regulation are charac- 
teristic of the so-called "late" LPS response [6-9]. 

Comparison of DC subset-specific LPS responses with 
published analyses of unfractionated DC responses 

In contrast to the extensive published microarray characterisa- 
tion of steady state DC subsets, most studies assessing LPS 
responses in DCs have used in vitro-derived DCs. We compared 
our data with these LPS response datasets to test how many of the 
significantly differentially expressed genes specific for each subset 
had also been identified within the published datasets. DC 
responses to LPS have usually been modelled using murine 
bone-marrow (BM)-DCs derived from in vitro stimulation of BM 
cells matured into DCs in 5-8 day cultures containing granulo- 
cyte-macrophage colony-stimulating factor (GM-CSF), alone or in 
combination with IL-4, IL-3, IL-6, and/or stem cell factor (SCF) 
[36,37]. While BM-DCs can be subdivided into CD24- and 
GDI lb-expressing subpopulations [36-38], they have been 
analysed by microarray only as a 'mixed' population. LPS 
responses are dynamically regulated over time and influenced by 
multiple factors including the type and dose of LPS [13,31,39]. We 
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Figure 2. Subset-specific LPS-induced gene signatures. (A-B) Dendritic cells were stimulated in vivo with an intraperitoneal injection of 25ug 
LPS and their activation confirmed by measuring upregulation of CD86 24 hours later. (C-E) LPS stimulated DC subsets isolated from a second cohort 
of mice (n = 1 0, spleen cells pooled) were magnetic bead enriched, sorted, analysed by RNA-Seq, and then compared to steady-state controls (n = 5, 
spleen cells pooled). (C) Differences in LPS-induced gene expression were visualised in a Venn diagram. (D-E) Heatmaps showing the expression of 
key immune effector genes (D) uniquely regulated in either CDS (left) or CDIIb DCs (right) or (E) similarly regulated in both subsets. Data are 
presented as fold changes (4-LPS/-LPS). 
doi:1 0.1 371 /journal.pone.01 0061 3.g002 



reanalysed 5 published microarray datasets, containing a total of 
10 timepoints (datasets 10-14 listed in Table SI, [10-14]), and 
identified 12,886 LPS responsive genes (p-value <0.05) in BM- 
DCs. 

We then performed 2 different analyses comparing our 24 hour 
timepoint data from sorted DC subsets with the published datasets. 
In the first, we included all 12,886 published LPS-responsive genes 
in the comparison. Of the 705 LPS responsive genes identified in 
our in vivo studies, 484 (69%) were also identified in at least one of 
the in vitro BM-DC datasets (Figure 3A). The high degree of 
overlap between our in vivo activated DCs and in vitro stimulated 
BM-DCs further validates the quality of our results and suggests 
that a major part of the observed transcriptional response in our ex 
vivo DCs was directly TLR4-mediated. However, 289 (60%) of 



these 484 commonly identified genes were regulated in a subset- 
specific manner in either CD8 (143 genes) or CDl lb (146 genes) 
DCs after in vivo LPS stimulation (Figure 3A). Next we repeated 
our comparison using only the 24 hour timepoint from the 
GSE17721 dataset [13], the sample that was the most consistent 
with our experimental setup. A similar trend was observed in this 
analysis, with 98 (52%) of the 189 genes identified in both studies 
being regulated in a subset-specific manner in our dataset 
(Figure 3B). 

We also reanalysed the only published microarray dataset 
derived from splenic DCs isolated directly from LPS-treated mice 
(dataset 15 listed in Table SI, [15]). These cells were stimulated 
for 6 hours in vivo with a combination of LPS and anti-CD40, and 
isolated on the basis of CD 1 1 c expression without fractionation 
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Table 1. GO terms annotating LPS-induced genes in CDS versus CDllb DCs. 
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GO term over-representation analysis of LPS-induced genes in CDS and CDllb DCs. 

1 . The ratio of odds (Odds-Ratio) that a GO term is enriched in the selected DC subset was calculated as the odds of a differentially expressed gene divided by the odds 
of a non-differentially expressed gene occurring in the GO term. 

2. P-values are adjusted to control for multiple comparisons. 
*denotes not significant (p>0.05). 

doi:l 0.1 371 /journal.pone.Ol 0061 3.t001 



into subsets [15]. Despite difTerences in stimulus and timepoint, 
190 of our 705 genes were also identified in this dataset 
(Figure 3C). Consistent with our observed subset specificity, 47% 
of the 190 were differentially regulated between subsets in our 
RNA-Seq analysis (Figure 3C). 

These results indicate a degree of DC subset-specificity within 
the LPS response that had not previously been apparent from 
analysis of unfractionated populations such as BM-DCs and 
CDl Ic-expressing splenic DCs. Many of the previously identified 
'LPS-responsive genes' may be regulated within only a subset of 
the total population, and novel subset-specific genes may be 
missed in such analyses. 

Splenic DC subsets share a common set of LPS- 
responsive transcription factors and signalling molecules 

Subset-specific transcriptional responses to LPS have been 
reported in non-DC cell types, but the mechanisms underlying this 
specificity have not been thoroughly addressed [40,41]. To identify 
molecular mechanisms underlying the observed transcriptional 
dififerences between the 2 closely related splenic DC subsets in our 
study, we employed a number of systems biology approaches. As a 
first step, we tested whether there were major primary signalling 
differences between the 2 subsets by looking for differential 
expression and activity of key LPS-response molecules, including 
transcription factors. 

We performed over-representation analysis on transcription 
factor-target LPS-induced genes in each subset, using a list of 



transcription factor-gene interactions based on experimental 
evidence and downloaded from innateDB [42] . Both subsets were 
significandy enriched for genes regulated by the transcription 
factors Cebpb, Irfl, IrfB, Jun, Nfkbl, Rela, and Spl, although the 
individual genesets were only partially overlapping between the 
subsets (Table 2). Genes regulated by the transcription factor Egrl 
were enriched in both subsets, although the effect narrowly failed 
to reach statistical significance in CD8 DCs (p-value 0.062). 
Importantiy, all these transcription factors were constitutively 
expressed by both DC subsets and are known to play key roles in 
LPS responses [6-9]. 

We also compared the relative expression of core LPS response 
molecules in the steady-state and after LPS stimulation. 48 core 
molecules were curated either as key canonical signalling 
molecules in LPS responses defined in multiple publications [6— 
9] and/ or identified within the KEGG and Reactome databases 
(Figwe S2). Importantly, these core molecules are also known to 
serve as essential mediators of many other immune signalling 
pathways and thus would be predicted to function as Hubs in both 
primary LPS-TLR4 and secondary signalling pathways. Of the 48 
core molecules, only 3 (Ticam2, Tlr4 and Ikbke) showed 
differential expression between the 2 subsets, and none were 
consistently differentially expressed both before and after LPS 
stimulation (Figure S2). The TLR4 signalling adaptor protein, 
Ticam2 (Tril) was significandy upregulated in CDllb DCs only 
before stimulation, whereas Tlr4 itself and the signalling molecule. 
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Figure 3. Population vs. subset-specific LPS-responses. (A) Data from published studies comparing in vitro stimulated GIVl-CSF derived BM- 
DCs +/- LPS stimulation (datasets 10-14 listed in Table SI) were reanalysed for differential gene expression. All LPS-induced genes identified in BIVI- 
DCs at any of the 10 timepoints (p-value <0.05) were compared to LPS-induced genes in our study and visualised in a Venn diagram. (B) Only genes 
differentially expressed in the 24 hour timepoint sample from the GSE17721 BiVl-DC dataset (dataset 13 listed in Table SI) were compared to LPS- 
induced genes in our study and visualised in a Venn diagram. (C) Genes identified as significantly differentially expressed (p-value <0.05) in splenic 
GDI 1c* DCs stimulated in vivo with LPS and anti-CD40 and isolated after 6 hours were compared to in vivo LPS-induced genes in our study (isolated 
after 24 hours) and visualised in a Venn diagram (dataset 15 listed in Table SI). 
doi:1 0.1 371 /journal.pone.01 0061 3.g003 



Ikbke, were significantly upregulated in the CD 1 1 b subset only 
after stimulation (Figure S2). 

In this analysis, clear subset-specific difTerences in key 
transcription factors and core signalling molecules could not be 
identified. Instead, these data support a model in which both DC 
subsets signal through a common set of molecules. 

DC subset-specific responses are 'fine-tuned' by distinct 
pathway modulators 

We next performed network analysis on the difiFerentiaUy 
expressed genes in each subset, in order to identify potential 
modulators of the subset-specific responses. As a first step, we 
uploaded our defined list of 48 core LPS response molecules into 
the immune database and analysis platform InnateDB, to generate 
a network containing these molecules and their first-order 
interacting partners. This identified 2279 interacting genes. We 



then filtered this network for interacting nodes (genes) that were 
significandy differentially expressed in the DC subsets. This 
filtered analysis identified multiple diflFerences in LPS-responsive 
genes that interact with, and are known to modulate, the function 
of the 48 core signalling molecules (Figure 4A, Table S4). CD8 
DCs uniquely regulated Anxa2, Atfi, Birc2, Cd81, Ctnndl, 
Ddx58, Dnajbl, Egrl, Fkbp5, Gadd45g, Hspalb, Ikzf4, 1112b, 
lUrl, Itgam, Jak3, Ksrl, Map2k6, Mef2c, Nfkbia, Peli2, Prdml, 
Pygl, Relb, Sertadl, Spib, Stat4, Tgm2, TnfaipS, Trafl, and 
Zfp36. CD lib DCs uniquely regulated Atxnl, Bcl211, Bcl2111, 
Cd209b, Cdknla, Ercl, Fosl2, Fthl, Gadd45a, Hmoxl, Idl, Id2, 
Igf2r, 1115, lUrap, Jdp2, Map4k4, Myold, Nfe2, Nlrpl2, Nos2, 
Notch 1, NotchS, Optn, Pld3, Plin2, Pnrcl, Prdx5, Rhbdf2, Rhoc, 
Sh3bp5, Slcl2a2, Snap23, Sod2, Tlr2, Tlr4, Tnf, Tribl, Trim29, 
Tuba8, and Usp2. 
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Table 2. Expression of canonical transcription factor-target TLR4-dependent pathways mediating the LPS response in CDS and 
CDIIb DCs. 



Transcription Factor 


CDS DCs 




CD11b DCs 






Pval^ 


Odds-Ratio^ 


Pval^ 


Odds-Ratio^ 


Cebpb 


3.6e-05 


6.2 


2.8e-06 


6.2 


Irfl 


4.6e-08 


7.4 


8.1e-12 


8.4 


IrfS 


4.3e-06 


2.6 


9.6e-13 


3.3 


Jun 


9.6e-06 


7.5 


2.7e-10 


11.2 


Nfkbl 


4.3e-06 


6.5 


1.2e-ll 


9.7 


Rela 


2.2e-1 1 


7.9 


9.6e-13 


7.3 


Spl 


9.6e-06 


8.2 


l.le-02 


3.5 


Egrl 


6.2e-02* 


1.3 


5.1e-03 


1.3 



Transcription factor-target over-representation analysis of LPS-induced genes In CD8 and CDllb DCs. 

1. The ratio of odds (Odds-Ratio) tfiat a transcription factor-associated patfiway is enriclied in the selected DC subset was calculated as the odds of differentially 
expressed genes being regulated by the transcription factor divided by the odds of non-differentially expressed genes being regulated by the same transcription factor. 

2. P-values are adjusted to control for multiple comparisons. 
*denotes not significant (p>0.05). 

doi:l 0.1 371 /journal.pone.Ol 0061 3.t002 



A limitation of this method for finding differentially regulated 
pathway modulators is that their identification was based on 
known interactions with a predefined list of core signalling 
molecules involved in LPS responses. As an alternate unbiased 
approach, we generated first-order interaction subset-specific 



networks from all differentially expressed genes in each subset, 
and then applied an unbiased subnetwork analysis. This success- 
fully yielded core subnetworks containing 3 1 7 individual nodes for 
CDS DCs and 297 nodes for CDl lb DCs (Figure 5, Figures S3- 
S4, Tables S5-S6). These unbiased subnetworks, which include 



Mouse Subset-Specific Expression of 
Tlr4 Pathway-Interacting Proteins 



° Thioglycolate Bone-Marrow 
CDS CDIIb IVIacrophages IVIacrophages 




Human Subset-Specific Expression of 
Tlr4 Pathway-Interacting Proteins 



V51 V62 L» Retinal Choroidal ■= Cord-Blood Cord-Blood 



y6T cells y6T cells Endothelial Cells Endothelial Cells Monocytes Neutrophils 




Figure 4. Subset-specific expression of pathway modulators. LPS-induced regulation (4-LPS/-LPS) of pathway nnodulators in nnouse (A-B) 
and hunnan (C-D) cell subsets, visualised in Venn diagrams. (A) CDS and CDl 1 b DCs from our study. (B) Thioglycolate-elicited peritoneal macrophages 
and bone-marrov*/ derived macrophages [26]. (C) V61 and V62 yS T cells [40]. (D) Retinal vascular endothelium and choroidal endothelial cells [41]. (E) 
Cord blood peripheral blood monocytes and neutrophils [48]. Datasets in (B-E) are listed as 16-19 in Table SI. 
doi:1 0.1 371/journal.pone.01 0061 3.g004 



PLOS ONE I www.plosone.org 



7 



June 2014 | Volume 9 | Issue 6 | e100613 



Subset-Specific Responses to LPS in DCs 



the most interconnected genes involved in the LPS response of 
each DC subset, demonstrated less than 50% overlap with each 
other, indicating significant differences in response between the 2 
subsets. To initially characterise the 2 subnetworks. Pathway (fi'om 
Reactome) and GO term over-representation analysis was 
performed on the genes comprising the subnetworks. No 
significant differences in Reactome Pathways were observed 
(Table 3). Both subnetworks were significantly enriched in nodes 
annotated by Reactome as being in pathways relating to TLR4 
signaUing, including both Myd88-dependent and -independent 
cascades (Table 3). GO terms relating to the LPS-response, the 
MAPK and NFkB cascades, and the general inflammatory 
responses were also over-represented in both networks (Table S2). 

Next we examined the extent of node interconnection in the 
unbiased subnetworks in order to identffy Hubs, which are defined 
as the most interconnected nodes within a network (nodes with the 
highest number of interactions). Hubs receive and integrate signals 
from multiple signal transduction pathways and are thus thought 
to be essential regulators or modulators of cell signalling. In the 
context of our in vivo LPS response model, the identified Hubs 
would be predicted to be integrating primar)- LPS-TLR4 and 
secondary immune signals to generate the observed immune 
outcomes. The degree of interconnection for each node was scored 
using the cytoscape plugin Cytohubba [43] . This analysis revealed 
that many Hubs (defined as nodes with more than 5 interactions) 
were the core LPS response molecules curated from the literature 
(Figure 5, Figure S2, Tables S5 and S6). Tlr4, Trafb, Ikbke, Irf7, 
Nfkbia, Fos,Jun, and Mapkl were Hubs in both CDS and GDI lb 
DC subnetworks (Figure 5). When overlayed on the linear TLR4 
KEGG pathway, these common Hubs mapped to both the 
Myd88-dependent and -independent pathways (Figure S5). This 
analysis further supports a model in which transcriptional 
responses to both primary and secondary signals are orchestrated 
via molecular pathways common to both subsets, consistent with 
the transcription factor and core signalling molecule expression 
data. 

However, we did identify some core signalling molecules that 
were selectively present in only one of the 2 unbiased subnetworks. 
Thus the CDS DC subnetwork uniquely contained Map3k7, 
Mapkl 1, and Tab2, while the GDI lb DC subnetwork uniquely 
contained Ikbkg, Iraki, MapkS, Myd88, Nflsbl, Rela, Ripkl, 
Tirap, ToUip, and TraB (Figure 5). Once again, these Hubs 
mapped to both the Myd88-dependent and -independent path- 
ways (Figure S5), potentially suggesting that modulation of a 
'common' signalling cascade occurs at different intervention points 
in the 2 subsets. 

Given the large number of predefined core LPS response Hubs 
identified within both unbiased subnetworks, we analysed the 
subnetworks (Figure 5, Tables S5-S6) for the presence of the 
pathway modulators identified in the initial filtered analysis 
(Figure 4A, Table S4). A high degree of overlap was observed 
between the two analyses. 83% of the GDI lb and 96% of the 
CD8 subset-regulated nodes identified in the filtered analysis were 
also identified in the unbiased subnetworks. In addition, many of 
these pathway modulators were among the most highly intercon- 
nected Hubs (defined as nodes with more than 15 interactions) 
within the unbiased subnetworks, or directiy interacted with the 
most highly interconnected Hubs, consistent with a key functional 
role. Atf3, Ep300, Gadd45g, Ikzff, Jak3, Relb, Tnfaip3 and 
ZbtblG, which were selectively present in the CD8 DC unbiased 
subnetwork, and Cdknla, Fthl, Gadd45a, Id2, Notch3, Ser- 
binb9/Spi6, and Tradd in the GDI lb DC subnetwork, are all 
known key modulators of cell signalling. Our analysis indicates for 



the first time that these molecules are also candidate modulators of 
DC subset-specific responses to LPS. 

We also identified a number of additional cell surface receptors 
including IL-7R and CXCR4 as Hubs (5-20 interactions) present 
within the subnetworks of both subsets (Figure 5). These receptors 
are known to recognise the LPS-inducible Ugands IL-7 and 
CXCL12, respectively [44-^7], and their presence within both 
signalling subnetworks is consistent with an important role for 
secondary immune mediators in influencing the observed 
transcriptional responses in the context of the immune microen- 
vironment. 

Based on these results, we suggest that LPS responses are 
regulated through a set of common pathway molecules in both DC 
subsets. Subset-specific responses are achieved by differential 
regulation of known pathway modulators that subsequently 'fine- 
tune' signalling by means of their interactions with common Hubs 
and multiple other signalling molecules involved in primary (LPS- 
TLR4) and secondary cytokine/chemokine/ceU-ceU interaction 
pathways (Figure 6). 

The model of subset-specific LPS response regulation is 
applicable to non-DC cell populations 

To test whether regulation of transcriptional responses to LPS 
in other cell subsets is consistent with our model, we reanalysed 
published datasets for which two distinct but related cell subsets 
were stimulated with LPS in the same experiment (datasets 16-19 
listed in Table SI). This meta-analysis included paired datasets 
comparing mouse thioglycate-elicited peritoneal macrophages 
with bone marrow-derived macrophages (Figure 4B, Table S7, 
[26]), human V51 versus V52 y5 T cells (Figure 4C, Table SB, 
[40]), human retinal vascular versus choroidal endothelial cells 
(Figure 4D, Table S9, [41]), and cord-blood monocytes compared 
with neutrophils (Figure 4E, Table SIO, [48]) These datasets all 
used in vitro stimulated cells and thus measured only primary 
TLR4-mediated outcomes. 

We repeated the filtered and unbiased network analyses, as 
described above, on datasets 16-19. Consistent with our DC data 
(Figure 4A), subset-specific gene regulation was seen in every 
paired dataset (Figure 4B-E, Tables S7-S10). Once again, there 
was a large overlap between the filtered and unbiased analyses, 
and many of the most interconnected Hubs in each of the 
unbiased subnetworks were either core TLR4 signalling molecules 
or candidate modulators identified in the original filtered analysis 
(Tables S11-S18). As predicted by our model, many of these 
modulators were selectively present in only one of two paired 
subnetworks (Tables SI 1-18). 

Discussion 

Dendritic cells are essential in triggering and tailoring the 
adaptive immune response. Understanding how defined DC 
subsets differentially respond to pathogenic signals is a crucial 
step in understanding how subset-specific control of immune 
outcomes is achieved. To address this question, we employed 
systems biology approaches to globally characterise the responses 
of splenic CD8 and GDI lb DC subsets following LPS stimulation. 
Our methodolog).- focused on identifying sets of modulators (the 
product of network analysis) rather than specific genes, and was 
specifically designed for th(; analysis of single pooled samples. As 
our conclusions are based on a self-reinforcing network analysis 
[3-5], combined with a hyper-stringent method for choosing 
differentially expressed genes [33], the need for multiple replicates 
was reduced compared to methodologies considering individual 
gene expression events. In addition, we used extensive meta- 
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Figure 5. Network analysis of LPS-responslve genes in CDS and CD11b DCs. Individual network analyses were carried out on the 
transcriptional response of (A) CDS and (B) CDIIb DCs stimulated in vivo with LPS as compared to steady-state. Subnetwork analysis was used to 
enrich networks in an unbiased manner for interactions with differentially expressed genes. The full subnetworks for each subset (Right Panels) were 
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edited to show a smaller core network for ease of visualisation while retaining the original topology. These smaller core networks include those genes 
with the largest Hub degree (interconnectivity with other genes in the full subnetwork), all core LPS response molecules, and key subset-specific 
modulators interacting with these genes. Subset-specific nodes are labelled in blue (unique to the CDS subnetwork) or orange (unique to the CD1 1 b 
subnetwork), while nodes labelled in black text are present in both subnetworks. The size of each node is proportional to its Hub degree, while node 
colour indicates relative gene expression (4-LPS/— LPS). Square nodes represent core LPS response molecules. The full subnetwork diagrams, made 
using the Cytoscape plugin Cerebral to show the cellular localisation of each gene, are displayed in Figures S3 and S4, while a list of nodes/network 
characteristics is provided in Tables S5 and S6. 
doi:1 0.1 371 /journal.pone.01 0061 3.g005 



analyses to validate the quality of our data in the absence of 
technical replication (Figures IC, 3, Table SI). 

Using RNA-Seq on cells isolated ex vivo from pooled biological 
replicate animals and immediately FACS sorted, we have shown 
that CD8 and GDI lb DC subsets respond uniquely following in 
vivo LPS stimulation (Figure 2). We have defined new subset- 
specific biomarkers and mechanistic information in addition to 
confirming previously identified differences in surface marker 
expression and cytokine secretion following LPS stimulation in 
vitro and in vivo [34,37,49]. One unique finding from our analysis 
was the unexpectedly wide range of cytokines/chemokines and co- 
stimulatory molecules differentially regulated in CD 1 1 b and CDS 
DCs (Figure 2C). This might indicate an inherent difiference in the 
use of soluble versus ceU-to-ceU dependent mechanisms for 
modulating immune responses in the 2 DC subsets. Our analysis 
also implicates differential responsiveness to external stimuli as one 
of the mechanisms by which DC subsets regulate their function- 
ality. To our knowledge, this is the first comprehensive examina- 
tion of the inflammatory transcriptional response in well-defined 
DC subsets in vivo. In contrast to a recent report in which steady- 
state plasmacytoid DCs, conventional DCs, and the CDS DC 
subset were individually subjected to network analysis of only those 
genes uniquely expressed by DCs, we included all expressed genes 



in a direct comparison of subsets, thus highlighting the behaviour 
of key functional molecules expressed in common with a range of 
other cell populations [20]. Several other groups have also 
identified core modulatory molecules within immune response 
pathways, but have focused on related cell types, rather than 
subsets within a single cell type [50-54]. 

A major implication of our study is that the transcriptional 
response to LPS is DC subset specific, and that analysis of 
unfractionated populations wUl by definition over-generalise such 
responses. To highlight this concept, we compared LPS responsive 
genes in our study to those identified in published datasets. As 
predicted, datasets from in vitro stimulated BM-DCs and in vivo 
stimulated unfractionated CDllc^ DCs included many LPS 
responsive genes that were regulated in a subset-specific manner 
within our study (Figirre 3A-C). This observation has major 
relevance to our understanding of responses to LPS, which has 
usually been defined by the study of mixed populations, and has 
often been generalised/extended across additional populations 
without experimental verification. Our study suggests that 
transcriptional control of signalling networks may need to be 
refined in a subset-specific context. 

While we observed a high degree of overlap with published 
datasets from LPS-stimulated BM-DCs, over 30% (221 out of 705) 



Table 3. Reactome pathways enriched within the CDS and CDllb DC subnetworks. 



Reactome Pathway 


CDS DCs 


CD11b DCs 




Pval 


Pval 


Activated TLR4 signalling 


1 .2e-06 


7.3e-l 1 


Adaptive Immune System 


1 .7e-04 


3.1e-09 


Cell-Cell communication 


2.5e-05 


4.4e-03 


Cytokine signaling in Immune system 


5.16-21 


2.4e-19 


Growth hormone receptor signaling 


3.6e-10 


1.6e-10 


Immune System 


2.9e-26 


4.6e-29 


Innate Immune System 


1.5e-11 


8.0e-16 


Integrin cell surface interactions 


3.0e-09 


1 .4e-08 


Interferon signaling 


8.3e-11 


3.4e-05 


lnterleukin-2 signaling 


1 .8e-06 


7.3e-06 


lnterleukin-3, 5 and GM-CSF signaling 


2.6e-07 


8.9e-06 


MyD88 dependent cascade initiated on endosome 


1 .4e-04 


2.0e-07 


MyD88-lndependent cascade 


2.3e-06 


1.2e-11 


NFkB and MAP kinases activation mediated by TLR4 signalling repertoire 


5.9e-05 


6.0e-09 


RIG-I/MDA5 mediated induction of IFN-alpha/beta pathways 


2.6e-07 


1.1e-07 


Signaling by Interleukins 


8.5e-13 


6.8e-18 


Toll Like Receptor 4 (TLR4) Cascade 


1 .2e-06 


9.6e-12 


Toll Receptor Cascades 


2.6e-07 


2.3e-13 


TRIF mediated TLR3 signaling 


6.0e-05 


6.0e-09 



Reactome pathway over-representation analysis of nodes within the CD8 or CD1 lb subnetworks. P-values are adjusted to control for multiple comparisons. 
doi:1 0.1 371/journal.pone.01 0061 3.t003 
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CDS DCs 



TLR4 



CD11b DCs 




Gene Expression 

Figure 6. The pathway modulation model. Schematic represen- 
tation of our proposed model. LPS-induced signaling occurs through a 
common set of signaling molecules. LPS induces the expression of 
subset-specific TLR4 pathway modulators, which 'fine tune' signaling 
and allow for distinct immune outcomes in closely related cell subsets. 
These pathway modulators likely integrate both signals derived directly 
from TLR4 and exogenous signals from the microenviroment that 
contribute to the subset-specificity of the response. 
doi:1 0.1 371 /journal.pone.01 0061 3.g006 

of the LPS responsive genes identified in our study liad not 
previously been identified in any BM-DC dataset (Figure 3A), 
while 73% were not identified in the BM-DC dataset from the 
same timepoint (24hr) as our study (Figure 3B). This may be due to 
many factors, including the influence of additional signals in the 
splenic microenvironment, differences in the starting DC popu- 
lation, experimental variables such as the use of microarrays and/ 
or differences in experimental set-up/reagents, or a combination 
of these factors. GM-CSF- induced BM-DCs are known to 
comprise a mixed population of 'inflammatory' DCs (CD24- 
and GDI lb-expressing subsets) that poorly reflect the in vivo 
steady-state DC populations in terms of surface marker expression 
and cytokine secretion [36-38]. Thus the comparison between 
BM-DCs and steady-state splenic DC subsets is complicated both 
by the poor correlation between populations and by the 'mixed' 
nature of the BM-DC population. Despite this, our results suggest 
the use of murine BM-DCs as model 'DCs' is failing to capture the 
subset-specificity of LPS responses within the DC compartment. 
Fms-like tyrosine kinase 3 ligand (Flt3L) may support a more 
physiological in vitro model for investigating DC signalling in a 
subset-specific manner, since Flt3L-generated BM-DCs more 
closely resemble the splenic CDS DC [36]. Similar to the BM- 
DC comparisons, the disparity between our findings and the 
published analysis of in vivo LPS-stimulated CDllc^ DCs 
(Figure 3C) may in part be explained by experimental differences 
including the addition of anti-CD40 to the LPS stimulus, 
combined with the much shorter stimulus time (6 hr versus 
24 hr) [15]. 

To explore how subset-specific differential responses to LPS 
stimulation may arise, we performed network analysis on LPS- 
regulated genes in the 2 DC subsets. We initially defined a set of 48 
core LPS response molecules based on their well established roles 
in propagating and modulating LPS responses. The lack of gross 
signalling differences between the subsets was indicated by 
comparable expression of the vast majority of the 48 core 
molecules in the 2 subsets (Figure S2). Unbiased subnetwork 
analysis revealed that many of the core LPS response molecules 
were also present as core Hubs within the unbiased subnetworks 
(Figure 5), while Reactome pathway analysis indicated that both 



the Myd88-dependent and -independent pathways were highly 
enriched within the both subnetworks (Table 2). Thus it appeared 
that primary LPS-dependent signalling was likely to be similar in 
the 2 DC subsets. One notable exception was the increased steady- 
state expression of the signalling adaptor Ticam2 (Tril) in CDl lb 
compared with CDS DCs (Figure S2), which might suggest a bias 
in CDl lb DCs towards signalling via the 'Myd88-independent 
pathway' mediated by Trif However overlaying the core LPS 
signalling molecules within the CD 1 1 b subnetwork (Ikbkg, Irak 1 , 
Mapk8, Myd88, Nflcbl, Rela, Ripkl, Tirap, Traf3, and ToUip) 
onto the KEGG TLR4 pathway (Figure S5) indicated a distribu- 
tion across both Myd88-dependent and -independent arms. This 
was also the case for the subset of core molecules common to both 
networks (Ikbke, Irf7, Fos, Jun, Mapkl, Nfkbia, and TralB), once 
again supporting a model in which both subsets share the same 
primary TLR4-dependent signalling machinery. The three core 
molecules uniquely present in the GD8 subnetwork (Map3k7, 
Mapkl 1, and Tab2) mapped to the Myd88-dependent pathway. 
However CD8 DCs have previously been shown to signal through 
Trif after polyLC stimulation, indicating that they also possess a 
functional Myd88-independent pathway [55]. 

The transcription factor-dependence of the differentially 
expressed genes was also comparable between the 2 DC subsets 
(Table 2). The identified transcription factors, including Nlkbl, 
Rela and Jun (AP-1), are all known to adopt key roles in functional 
aspects of the LPS response [6-9]. Although Irf8 was expressed at 
consistently higher levels in GD8 DCs (Table 2, Figure IB) and is 
known to be critical for CD8 but not CD lib DC development 
[56], genes regulated by IrfB were highly enriched in both CDl lb 
and CD8 subsets, consistent with the known impairment of 
responses to CpG and LPS in IrfS-deficient CDl lb DCs [56]. 
Collectively, these data suggests that both DC subsets utilise a 
common core pathway, characterised by a well established set of 
signalling Hubs and transcription factors. 

The absence of clear differences in known pathway components 
supports a more subtie mechanism of signal modulation (Figure 6). 
Both filtered and unbiased network analyses clearly identified sets 
of genes with known regulatory or modulatory function uniquely 
present within the subset-specific subnetworks (Figures 4A, 5, 
Tables S4-S6). We hypothesize that LPS-dependent transcrip- 
tional regulation of these genes mediates the observed subset- 
specific differences in the response to the combination of primary 
(LPS-TLR4) and secondary immune signals to which the DC 
subsets are exposed in vivo. Within CD8 DCs, the GO term 
'negative regulation of the inflammatory response' was over- 
represented, supported by the identification of known negative 
regulators AtB, Tnfaip3 (A20), and Zbtbl6 (PLZF) as key hubs in 
the network analysis, and suggesting that negative regulation may 
be a hallmark of the CD8 subset [57-59]. Interestingly, Atf3 is a 
predominant negative regulator of inflammation [51] and has 
been previously identified as one of the most important molecules 
identified within the signalling network of LPS-stimulated bone- 
marrow macrophages [54]. Differential regulation of MAPK 
signalling might by another mechanism by which DC subsets 
mediate subset-specific responses, as Gadd45g (uniquely identified 
in the CD8 subnetwork) and Gadd45a (uniquely identified in the 
CDl lb subnetwork) are known to differentially modulate MAPK 
pathways in hepatoma cells [60]. Similarly, the CDl lb subnet- 
work uniquely contained TRADD, which is known to play a 
differential role in regulating NFkB and MAPK signalling in 
fibroblasts versus macrophages and is a likely key mediator of the 
unique CDl lb DC LPS-response [61]. Cdknla (p21/WAFl/ 
CIPl), also exclusively present within the CDl lb subnetwork, was 
another highly interconnected Hub that is known to be essential 
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for regulating LPS activation in macrophages [62,63]. Ttius our 
analysis has revealed previously unappreciated subset-specific roles 
for several known signaUing modulators. 

We propo.se that our model in which the fme-tuning of central 
immune pathways mediates subset specific responses is indeed 
relevant in a physiological context since our data were obtained 
from cells taken directly from mice with minimal manipulation. 

Our experimental design characterising 'late' in situ DC subset 
responses to LPS provides a more physiological model than an ex 
vivo stimulation assays. However, while both DC subsets can 
respond directly to LPS [31], our analysis is complicated by the 
integrated response of the DC subsets to additional exogenous 
signals arising from the splenic microenviroment. In agreement 
with this, we identified differential expression of a number of 
additional cell surface receptors within the unbiased subnetworks 
of the 2 DC subsets (Figure 5). This suggests that both DC subsets 
are integrating direct (LPS-TLR4) signals with those from other 
secondary immune mediators such as IL7 (via IL7R) and 
CXCL12 (via CXCR4), which are known to play multifunctional 
roles in regulating DC function [44—47]. However, the greatest 
strength of unbiased subnetwork analysis is the ability to identify 
key functional molecules (Hubs) within complex signalling 
networks, without relying on previously defined linear pathways. 
Therefore, the core LPS response molecules and pathway 
modulators identified here as Hubs are likely responsible for 
integrating signals derived direcdy from TLR4 with additional 
exogenous signals, resulting in the observed subset-specific 
responses. Indeed a major strength of this study is the physiolog- 
ical, complex nature of the in vivo stimulus, and its potential to 
serve as a basis for identifying key secondary immune signals 
regulating DC function following LPS exposure (such as CXCL12- 
CXCR4). However, further studies are needed to fully characterise 
how the observed subset-specific responses are controlled in the 
context of the immune microenvironment. 

While we have identified many candidate genes potentially 
involved in regulating subset-specific responses (as discussed 
above), our hypothesised model is primarily based on sets of 
subset-specific modulators regulating common signalling path- 
ways, rather than these individual genes. This ensures that our 
conclusions are relatively resistant to errors introduced by using a 
single pooled sample for each condition. As with any transcrip- 
tional-based study, complex functional studies, such as those using 
knockout/knockin models, will be required to validate individual 
genes that are actively regulating subset-specific function. How- 
ever, we believe that our approach represents an important first 
step in understanding how cell subsets may regulate their 
responses to common stimuli. 

Given our focus on sets of modulators, we chose to validate the 
suitability of our hypothesised model of pathway modulation and 
signal integration, rather than the expression of individual genes. 
To do this, we tested how well our model fitted to published 
datasets for both mouse and human cell subset responses to LPS. 
In each dataset, we identified cell subset-specific modulators that 
are known to interact with core LPS response molecules and 
which uniquely act as Hubs within their respective signalling 
networks (Tables S7-S10). While we cannot fully exclude other 
mechanisms contributing to the observed subset differences in 
these published datasets, the results strongly support our current 
model. Thus pathway modulation appears to represent a global 
mechanism allowing for tighdy controlled and specific responses to 
LPS in related but distinct cell populations. 

While we have subdivided DCs in this study on the basis of CDS 
and CD lib expression, multiple reports have suggested that the 
splenic DC network is much more complex and that splenic DCs 



can be further subdivided based on the expression of many 
additional surface markers including ESAMl [64], DCAL2 [65], 
CD207/Langerin [66], CD103 [67] and/or CD205/DEC-205 

[68]. Our model would predict that further subdivision on the 
basis of these markers would reveal additional layers of complexity 
in the subset-specific regulation of LPS responses. Further 
complexity also arises from the dynamic nature of LPS responses 
over time [13,39,53]. Previous network analysis of time course 
responses has shown that the majority of Hubs are only transiently 
involved as key regulators under certain conditions, and that even 
permanent Hubs redefine their interactions dynamically [5] . Since 
we consistentiy identified differential subset-specific regulation of 
TLR4-interacting proteins in datasets from multiple timepoints 
and cell populations in mouse and human (Figure 4, Tables S7— 
S 1 0), it is likely that this means of signalling modulation and fine 
tuning plays a critical role throughout the response. 

Materials and Methods 

Mice and treatment 

All mice were housed under specific pathogen-free conditions in 
the Centenary Institute Animal Facility. [C57BL/6x B10.BR]F1 
mice on a CD45.1/CD45.2 heterozygous background were used 
for all experiments. Our unpublished results have identified no 
differences in phenotype and function of splenic DC subsets 
isolated from [C57BL/6x BlO.BRjFl mice and their C57BL/6 
counterparts. Mice were injected intraperitoneally with 25 ng LPS 
(rough strains from Salmonella enterica serotype Minnesota Re 595, 
Sigma-Aldrich) per mouse. Animals were sacrificed mice 24 hours 
after LPS injection. 

Ethics Statement 

Approval for all animal experimentation was obtained from the 
Animal Ethics Committee at the University of Sydney. 

Flow cytometry and cell sorting 

Spleens were pooled from 5 (control) or 10 (LPS-stimulated) 
animals, digested with 2 mg/ml CoUagenase IV from Clostridium 
Mstxtlytkum (Sigma-Aldrich) and a single cell suspension prepared as 
described previously [69]. Cells for RNA-Seq were prepared by 
cell sorting. Briefly, splenocytes were stained for B220 (clone RA3- 
6B2) prior to DC selection. DCs were isolated with anti-CD 11c 
MicroBeads (clone N418, MUtenyi Bio tec) after enrichment for 
Lineage (Terl 19/B220/CD19/CD3/Gr-l)-negative cells using 
rat anti-mouse antibodies and anti-rat IgG MicroBeads (Miltenyi 
Biotec). Selected DCs (78% purity) were stained for GDI Ic (clone 
HL3), GDI lb (clone Ml/70), CDS (clone 53-6.7) (all from BD/ 
Pharmingen) and pan-MHCII (clone M5/114, eBioscience). 
Staining was performed in PBS containing 5% ECS and lOmM 
EDTA, non-specific staining was blocked with unconjugated anti- 
CD16.32 (clone 2.4G2) and DAPI was used to exclude dead cells. 
B220"MHGII+CD1 lc+ cells were sorted into GDI lb+ and CDS"^ 
subsets on Aria-Ilu (BD) to 99.1'/'o and 86.3% purity, respectively). 
For analysis, DAPI-negative (live cells) events were gated on 
forward scatter height vs. area to exclude doublets. CD86 
expression was detected using anti-CD86 (clone GLl, BD/ 
Pharmingen). To identify monocyte-derived DCs, cells were 
additionally stained with anti-FcsRla (clone MAR-1, 
eBioscience), GD64/FcyR1 (clone X54-5/7.1, BD/Pharmingen) 
and Gr-1 (unconjugated, detected with anti-Rat- Alexa488, Invi- 
trogen). 
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RNA isolation and RNA-Seq 

Total RNA was isolated from sorted DCs (>98% purity in 
every sample) using RNeasy Micro Kits (Qiagen). A minimum of 
200,000 cells and 100 ng of RNA per condition were used for 
RNA-Seq. RNA quality was determined using a Bioanalyzer 2100 
(Agilent Technologies), and cDNA libraries were prepared from 
total RNA according to the lUumina TruSeq RNA sample 
preparation guide. Unique adapter indexes (lUumina) were 
attached during sample pr(;p and samples were run pooled and 
loaded into a single flow cell lane to reduce technical variability. 
RNA-Seq was performed on a GAIIx instrument (lUumina), using 
a single read run with 36 amplification cycles (29 sample -h7 
adapter/index sequence). 

Data Processing 

Raw basecall data was converted to FASTQ^ sequence files 
using OfT-Line Basecaller (Ilumina) and a custom Perl script. 
Reads were aligned to the mm9 mouse genome with TopHat 
version 2.05 and Bowtiel version 0.12.7 [70]. Reads were initially 
mapped to Ensembl transcripts with the search for novel junctions 
disabled, using standard Tophat filtering/ stringency parameters 
[70]. Genomic coordinates were then transformed into counts of 
protein-coding Ensembl genes. To do this, a chimeric gene-model 
was first defined by merging all protein-coding transcripts for a 
given gene. Transcripts that had reads in less than 50% of their 
exons in all samples were defined as not expressed and were 
excluded from the chimeric transcriptome. Reads that overlapped 
the chimeric genes were counted using the htseq-count script in 
the intersection-nonempty mode (http://www-huber.embl.de/ 
users/anders/HTSeq/doc/count.html). The script discards 
multi-mapped reads as well as reads that overlap multiple distinct 
genes, to generate a file of uniquely mapped gene counts. No 
additional gene filtering was performed. Total mapped reads to 
protein coding genes were 6702672 for steady state and 5536485 
for LPS-stimulated CDS DCs, and 6368840 for steady state and 
2582238 for LPS-stimulated CDllb DCs. 

Transcriptional Analysis 

Analysis of RNA-Seq gene count data was performed in R using 
Bioconductor ([71] r02290). Differential gene expression was 
calculated using the DESeq package [33] . DESeq was run with the 
method = blind and sharingMode = fit-only settings for single- 
replicate experiments. This sets all samples as replicates of each 
other when calculating the variance. This method tends to be 
over-conservative as compared to a replicated experiment ([33] 
r02353). Genes with an associated p-value £0.05 were scored as 
differentially expressed between samples. A fuU list of expressed 
genes is provided in Table SI 9. Gene ontology (GO) and 
transcription factor-target over-representation analyses were cal- 
culated using the Wallenius distribution in the goseq package, 
which normalises for RNA-Seq length biases ([72] r02354). Tests 
in which fewer than 10 genes in a term were observed in both 
subsets were excluded from further analysis. The threshold of 
significance in the ORA tests was defined as a Benjamini- 
Hochberg [73] adjusted p-value £0.05. Odds-Ratios were 
calculated as the odds of a differentially expressed gene occurring 
in the ORA category divided by the odds of a non-diflferentially 
expressed gene occurring in the ORA category, so that numbers 
greater than 1 are considered to reflect associations that are likely 
to be real. GO-Terms to gene mappings and Reactome Pathway 
to gene mappings were obtained using the biomaRt package [74] . 
Transcription factor to gene mappings were downloaded from 
InnateDB after searching for protein-gene interactions between all 
genes (including those predicted by orthology) [42] . RNA-Seq data 



was deposited in the NCBI Gene Expression Omnibus (GEO) 
repository (GSE42573). 

Network analysis 

DiflFerentially expressed (DE) genes from multiple analyses were 
uploaded separately into InnateDB [42], a specialised interactome 
database containing all known protein-protein interactions but 
highly curated for immune protein reactions, to generate a list of 
interactions between DE genes in the dataset, and with first-order 
non-DE proteins with curated experimental evidence of an 
interaction. Interactions predicted l)y orthology were included in 
all analyses. To analyse the complex functional relationships 
between genes comprising the observed LPS stimulation signa- 
tures, regulated genes were viewed as nodes in a network 
(■'network analysis"), connected to one-another by their protein- 
level interactions (edges), as previously described [75,76]. Briefly, 
interaction networks were visualised in Cytoscape [77] using the 
Cerebral plugin to show gene location relative to the ceU [78]. 
These networks were filtered for protein-protein interactions after 
removing duplicate edges, self-loops and the general ubiquitin, 
Ubc. Ubc interacts with ~3000 proteins and its inclusion thus 
biases subsequent subnetwork analyses. Subnetwork analysis was 
carried out on each of the networks using the jActive plugin for 
Cytoscape, and the top significant subnetworks were ranked on the 
basis of their calculated Z-score [79,80]. Multiple significant 
subnetworks were merged for some comparisons based on their 
high degree of overlap and associated z scores within a stratum. 
Ilighl)' interconnected gene nodes within these subnetworks are 
referred to as Hubs, and represent key molecules involved in signal 
trafficking. While the term "Hub" cannot be rigorously defined in 
the context of computational network analysis, we have defined 
Hubs as nodes with 5 or more interactions and the "most highly 
interconnected" Hubs as nodes with 15 or more interactions. Hub 
degree was scored using the Cytoscape software plugin cytoHubba 
with the higher scoring nodes predicted to represent essential 
molecules [43]. 

In separate filtered analyses, core signalling molecules in the 
TLR4 pathway in mouse and human were uploaded into 
InnateDB to generate a list of 2279 first-order interacting proteins 

[42] . This network was filtered for interacting nodes significantiy 
differentially expressed in a cell subset and visualised as a venn 
diagram representing each interacting nodes' differential expres- 
sion in one or both of the subsets. 

Meta-Analyses 

Normalised datasets were downloaded from NCBI GEO using 
the Bioconductor package GEOquery [81]. GSE15907 and 
GSE32381 were downloaded manually from NCBI GEO as raw 
CEL files, quantile normalised and RMA background corrected. 
Datasets were included in our analyses if they contained at least 
two biological replicates (see Table SI for a list of reanalysed 
datsets). Differential expression was calculated using the limma 
package. Genes were defined as significant at a Benjamini- 
Hochberg [73] adjusted p-value <0.05. A minimum fold change 
cut-off of 2 was applied for network analysis. Genes from yS T cells 
(GSE3720), endothelial ceUs (GSE7850), and cord blood mono- 
cytes and neutrophils (GSE39840) were defined significant at a 
non-adjusted p-value cut-o£Fof 0.05 and a fold-change cutoff of 2. 
To compare our RNA-Seq study to previously published 
microarray datasets characterising steady-state splenic DC subsets 
(datasets 1-9 listed in Table SI), we used a hypergeometric test 
approach. To do this, we first identified those genes in each dataset 
that were significantly difierentially expressed in CDS versus 
CDllb DCs at a p-value cut-ofi" of 0.05. We then used a 
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hypergeometric test to calculate the significance of overlap 
between each of these gene lists and the gene list derived from 
our RNA-Seq study. 

Supporting Information 

Figure SI DC purification strategy excludes monocytes 
and inflammatory monocyte-derived DCs. Spleen cells 
from control and LPS-injected mice were subjected to the pre-sort 
bead selection procedure as described in Materials and Methods. 

Selected cells were stained for MHCII, CDllc, CD64/F(:yRl, 
FcsRla, Gr-1 and GDI lb, and analysed for the presence of 
contaminating monocytes, inflammatory monocytes and mono- 
cye-derived DCs. (A) The lineage (CD 19, B220, CDS, Gr-1, 
TerlI9)-negative MHCII''"CD1 Ic"*" gating strategy excludes 
FcyRl^ceRIa^ monocyte-derived DCs. (B) Conversely, mono- 
cytes and inflammatory monocytes expressing Gr-1 do not 
significandy contaminate the MHCII^CD 1 1 c^ sorting gate shown 
in the right panels. 
(TIFF) 

Figure S2 Expression of core LPS response molecules. 

Comparison of core LPS response molecules in GD8 and CD lib 
DCs in the steady-state (—LPS) and after LPS stimulation (+LPS). 
Data are presented as fold changes (GDI lb/CD8). * Significandy 
differentially expressed before LPS stimulation; ** Significandy 
differentially expressed after LPS stimulation. 
(TIFF) 

Figure S3 Network analysis of LPS-responsive genes in 
CDS DCs. A network analysis was carried out on the 
transcriptional response of CDS DCs stimulated in vivo with 
LPS as compared to steady-state. Subnetwork analysis was used to 
enrich networks in an unbiased manner for interactions with 
differentially expressed genes. The figure was made using the 
Cytoscape plugin Cerebral to show the cellular localisation of each 
gene. The size of each nod(' is proportional to its Hub degree 
(interconnectivity with other genes), while node colour indicates 
relative gene expression (+LPS/ — LPS). Square nodes represent 
core LPS response molecules. Nodes labelled in blue text are 
present in the CD8 but not GDI lb DC subnetwork, while nodes 
labelled in black text are present in both. Networks were organised 
using the Cytoscape plugin Cerebral, which organises nodes based 
on their relative cellular location. For visualisation, only selected 
nodes are labelled. The fuU list of nodes/network characteristics is 
provided in Table S5. 
(TIFF) 

Figure S4 Network analysis of LPS-responsive genes in 
CDllb DCs. A network analysis was carried out on the 
transcriptional response of GDI lb DCs stimulated in vivo with 
LPS as compared to steady-state. Subnetwork analysis was used to 
enrich networks in an unbiased manner for interactions with 
differentially expressed genes. The figure was made using the 
Cytoscape plugin Cerebral to show the cellular localisation of each 
gene. Node size is proportional to its Hub degree (interconnectivity 
with other genes/nodes), and node colour indicates relative gene 
expression (+LPS/ — LPS). Square nodes represent core LPS 
response molecules. Nodes labelled in orange text are present in 
the GDI lb but not GD8 DC subnetwork, while nodes labelled in 
black text are present in both. Networks were organised using the 
Cytoscape plugin Cerebral, which organises nodes based on their 
relative cellular location. For visualisation, only selected nodes are 
labelled. The full list of nodes/network characteristics is provided 
in Table S6. 
(TIFF) 



Figure S5 Subset-specific Hubs in relation to a KEGG 
pathway map of TLR signalling. Gore LPS response Hubs 
identified in the subnetwork analysis of CD8 or CDllb are 

identified by coloured dots and gene names (italics) overlayed on a 
KEGG pathway map. Black dots and text indicate Hubs identified 
in both subnetworks, blue indicates Hubs identified only in the 
CD8 subnetwork and orange indicates Hubs identified only in the 
GDI lb subnetwork. 
(TIFF) 

Table SI List of reanalysed datasets and their associ- 
ated references. 

(DOCX) 

Table S2 GO term over-representation analysis on 
nodes within the CDS or CDllb subnetworks. P-values 
are adjusted to control for multiple comparisons. 

(CSV) 

Table S3 Differentially-expressed genes identified by 
comparing LPS stimulated with steady state expression 
data for each DC subset. 

(XLSX) 

Table S4 Gene list of differential pathway modulators 
in CDS and CDllb DCs from this RNA-Seq study, as 
depicted in Figure 4A. 

PCLSX) 

Table S5 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in CDS DCs. 

(XLSX) 

Table S6 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in CDllb DCs. 

PCLSX) 

Table S7 Gene list of differential pathway modulators 
in thioglycolate-elicited peritoneal macrophages and 
bone-marrow derived macrophages, as depicted in 
Figure 4B. 

(XLSX) 

Table S8 Gene list of differential pathway modulators 
in V61 and V82 78 T cells, as depicted in Figure 4C. 

(XLSX) 

Table S9 Gene list of differential pathway modulators 
in retinal vascular endothelium and choroidal endothe- 
lial cells, as depicted in Figure 4D. 

PCLSX) 

Table SIO Gene list of differential pathway modulators 
in cord blood monocytes and neutrophils, as depicted in 
Figure 4E. 

(XLSX) 

Table Sll Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in thioglycolate-elicited macrophages. 

(XLSX) 

Table S12 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in bone-marrow derived macrophages. 

PCLSX) 
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Table S13 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in V81 T cells. 

(XLSX) 

Table S14 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in V82 78 T cells. 

(XLSX) 

Table S15 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in retinal vascular endothelial cells. 

PCLSX) 

Table S16 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in choroidal endothelial cells. 

(XLSX) 

Table S17 Full node lists and corresponding network 
characteristics for the subnetwork of LPS-responsive 
genes in cord blood monocytes. 

PCLSX) 
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